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L  SUMMARY  1  | 

This  final  report  summarizes  the  achievements  and  activities  on  the  AFOSR  Grant 
No.  83-0071  from  February  15,  1983  to  August  15,  1988. 

The  research  under  this  grant  has,  to  date,  resulted  in  four  journal  publications  and 
seven  conference  papers.  One  additional  paper  has  been  accepted  for  presentation  at  the 
AIAA  27th  Aerospace  Sciences  Meeting,  and  another  one  will  be  submitted  for  presen¬ 
tation  at  the  next  AIAA  Fluid  Dynamics  or  CFD  conference.  It  is  anticipated  that  there 
will  be  two  more  journal  publications  for  all  the  worit  that  has  been  completed.  Nine 
students  have  been  supported  fully  or  partially  on  this  grant  Two  graduated  with  a 
PhJ3.,  four  with  an  M.S.,  and  the  other  three  continue  to  work  on  their  dissertation 
topics. 

The  main  theme,  as  originally  proposed,  has  been  adaptive  solution  refinement  A 
novel  approach  called  Truncation  Error  bijection  (TEI)  was  introduced  during  the  course 
of  research.  The  idea  behind  TEI  is  very  simple,  i.e.,  the  exact  nodal  value  of  the  solu¬ 
tion  to  a  differential  equation  could  be  obtained  on  any  grid  and  from  solving  a 
difference  equation  that  models  the  differential  equation  if  the  truncation  error  were 
known.  Although  the  TE  is  not  known  in  general,  it  can  be  approximated  on  a  local  grid 
patch.  This  approach  of  approximating  the  local  error  due  to  discretization  in  effect 
decouples  a  problem  of  multiple  disparate  length  scales  into  problems  of  single  length 
scale  so  that  they  can  be  solved  more  efficiently  on  a  computer  than  the  original  prob¬ 
lem.  Three  types  of  applications  have  been  demonstrated.  In  addition  to  solution 
refinement  by  TEI,  we  have  shown  that  the  decoupling  of  the  unsteady  computation 
from  the  steady  one  by  TEI  could  significantly  reduce  the  computing  time  and  storage 
for  flutter  prediction,  and  that  viscous  effects  can  be  computed  separately  and  injected 
into  the  solution  of  an  inviscid  solver  for  viscous  flow  computation.  Some  of  the  advan¬ 
tages  of  this  approach  are:  it  requires  very  little  modification  to  the  base  solver;  no  com¬ 
patibility  problems  in  using  different  grids  and  different  solvers;  readily  suited  for 
multi-processors. 

This  method  has  also  been  applied  to  study  problems  related  to  the  dynamic  stal¬ 
ling  of  an  airfoil.  A  key  question  in  dynamic  stall  is  when  and  where  does  the  boundary 
layer  separate?  Our  method  enables  one  to  resolve  the  details  of  a  local  flow  field  such 


t 


as  the  behavior  of  the  boundary  layer  at  the  leading  edge  of  the  airfoil.  This  extension 
of  our  work  will  be  reported  in  detail  for  a  different  AFOSR  grant 

n  SYNOPSES  of  JOURNAL  PUBUCATIONS 

Four  journal  publications  are  attached  as  Appendix  1.  The  first  paper,  entitled 
"Adaptive  Refinement  with  Truncation  Error  Injection,"  outlines  the  TEI  methodology 
and  demonstrates,  with  model  problems  in  fluid  mechanics,  the  versatility,  efficiency  and 
accuracy  of  the  methodology.  The  second  paper,  entitled  "Refined  Numerical  Solution  of 
the  Transonic  Flow  Past  a  Wedge",  shows  an  application  of  TEI  on  one  of  the  classical 
problems  in  transonic  flows.  A  solution  accurate  to  1%  of  the  drag  value  of  the  flow 
over  a  wedge  has  been  obtained  using  TEI.  An  analysis  of  this  result  confirmed  for  the 
first  time  the  existence  of  a  weak  oblique  shock  right  after  the  expansion  at  the  shoulder 
in  addition  to  the  normal  shock  downstream.  The  third  paper,  entitled  "Computation  of 
Unsteady  Transottic  Aerodynaitucs  with  Truncation  Error  Injection,"  extends  the  TEI 
methodology  to  decouple  the  computations  of  the  unsteady  aerodynantics  due  to 
unsteady  body  motion  and  the  steady  aerodynantics  due  to  body  thickness  such  that  the 
unsteady  computation  can  be  done  more  efficiently  on  a  much  coarser  grid  than  that 
used  for  the  steady  flow  computation.  A  factor  of  sixteen  in  the  saving  of  computing 
time  has  been  demonstrated  using  examples  of  oscillating  airfoil  in  the  transonic  speed 
range.  The  fourth  paper,  entitled  "Computation  with  Error  Injection,"  was  an  invited 
paper  presented  at  the  Sixth  International  Conference  of  Numerical  Modelling  in  Sci¬ 
ence  and  Technology.  This  paper  reviews  the  TEI  methodology  and  further  generalizes 
it  to  accommodate  the  use  of  different  equations  on  different  grids. 

m  PAPERS  to  be  PUBUSHED 

A  paper  entitled  "Viscous-Inviscid  Interaction  and  Local  Grid  Refinement  Via 
Truncation  Error  Injection,"  has  been  accepted  for  presentation  at  the  AIAA  27th 
Aerospace  Sciences  Meeting  of  January,  1989  in  Reno,  Nevada.  It  will  be  shown  that 
accurate  prediction  of  the  flow  over  an  airfoil  can  be  obtained  by  solving  the  Euler 
equation  on  a  relatively  coarse  global  grid  with  viscous  effects  computed  separately  on  a 


boundary-layer  type  grid  and  injected  into  the  global  grid  solution  as  a  combination  of 
vorticiOr  and  truncation  error.  In  addition  to  the  efficiency  and  accuracy  gained  by  using 
TEI,  the  solution  on  the  local  grid  reveals  details  of  the  shock  structure  and  a  jet-like 
flow  emanating  from  the  root  of  the  normal  shock  in  the  shock  boundary  layer  interac¬ 
tion  zone.  This  result  has  already  shed  some  light  on  the  possible  mechanisms  causing 
the  onset  of  separation  of  dynamic  stall.  It  is  anticipated  that  a  further  extension  of  the 
TEI  methodology  to  equation  solvers  in  body-fixed  coordinates  and  nonstationary  grid 
will  allow  us  to  explore  this  important  area  of  research. 

Another  paper  under  preparation  for  presentation  and  publication  is  called  "An 
Efficient  Scheme  for  Three-Dimensional  Unsteady  Transonic  Computations".  Many 
techniques,  including  TEI,  have  been  incorporated  into  an  full  potential  code  in  general¬ 
ized  coordinates  to  achieve  both  accuracy  and  computational  efficiency  for  flutter  pred¬ 
iction  for  the  transonic  regime.  A  full  scale  application  of  this  scheme  on  the  prediction 
of  the  flutter  boundary  of  an  AGARD  standard  wing  will  be  demonstrated. 

These  papers  will  be  sent  to  AFOSR  as  soon  as  they  have  been  published. 

rv  THESIS  TITLES  and  AUTHORS 

"Numerical  Studies  of  Shock  Wave  Resolution  By  Mesh  Refinement",  Masters 

Report,  1984  —  J.  M.  Tripp 

"Computations  of  Unsteady  Transonic  Aerodynamics  with  Truncation  Error  Injec¬ 
tion”,  Masters  Report,  1985  ~  J.-K.  Fu 

"Refined  Numerical  Solutions  of  The  Transonic  Flow  Past  a  Wedge",  Ph.D. 

Dissertation,  1985  -  S.-M.  Liang 
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V  CONFERENCE  PRESENTATIONS 
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Liang),  Paper  No.  AIAA-85-1593,  presented  at  AIAA  18th  Fluid  Dynamics  and 
Plasma  Dynamics  and  Lasers  Conference,  Cincinnati,  Ohio,  July  1985. 
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COMPUTATION  WITH  ERROR  INJECTION 


K.-Y.  Fung  and  Brian  D.  Goble 


University  of  Arizona.  Aerospace  and  Mechanical  Engr.  Dept.,  Bldg.  16,  Tucson, 
AZ  83721  USA 


Abstract.  In  the  context  of  computational  fluid  dynamics,  the  adaptive  refinement 
technique  with  truncation  error  injection  due  to  Fung  et  al  ( 1984)  is  extended  to 
include  differences  in  the  governing  equations  being  used  on  the  fine  and  coarse 
grids.  This  method  allows  the  decoupling  of  a  problem  of  multiple  length  scales 
into  problems  of  simple  length  scale  which  can  be  more  efficiently  resolved  on 
separate  grids  and  using  different  formulations  pertinent  to  the  length  scales  than 
on  a  single  grid.  Three  variants  of  this  method  are  exemplified  with  applications  to 
problems  in  computational  fluid  dynamics.  Substantial  savings  in  computer  time 
and  storage  are  achieved  in  the  examples. 


Keywords.  Solution  refinement:  Truncation  error:  .Multiple  grid:  Multigrid. 


INTRODUCTION 

In  the  following  we  will  introduce  a  methodology 
which  allows  the  decoupling  of  a  complex  problem 
of  distinct  length  scales  into  problems  of  simple 
length  scale  so  that  they  can  be  solved  more 
efficiently  on  a  computer.  The  examples  we  have 
chosen  here  to  apply  this  methodology  to  are  prob¬ 
lems  in  computational  aerodynamics,  since  they  are 
of  our  prime  interest.  We  hope  a  review  of  our 
works  in  the  development  of  this  methodology  will 
make  it  known  to  other  disciplines  and  further  its 
applications. 

Computation  of  flows  over  aerodynamic  bodies  has 
reached  a  stage  where  solutions  to  most  engineering 
problems  can  be  found  with  some  degree  of  accu¬ 
racy.  Solvers  are  available  for  equation  sets  ranging 
from  the  Laplace  equation  for  subsonic  or  super¬ 
sonic  flow  to  the  Reynolds-averaged  compressible 
Navier-Stokes  equations  for  transonic  flow.  The 
choice  of  equations  to  be  solved  to  predict  some 
given  flow  situation  depends  on  such  variables  as 
degree  of  accuracy  required,  computer  resources 
available,  codes  available,  complexity  of  geometry, 
and  complexity  of  flow  situation.  Ideally,  the  most 
accurate  code  available  would  be  used  on  a  grid  with 
sufficient  resolution  to  capture  all  of  the  relevant 
physics  of  the  flow  field.  However,  the  computer 
resources  in  terms  of  memory  and  CPU  time 
required  to  implement  the  desired  code  on  a  suit¬ 
able  grid  are  always  limited. 

The  generation  of  such  a  grid  is  one  of  the  main 
stumbling  blocks  in  the  solution  process.  One  way 
to  simplify  the  grid  generation  process  and  to 
increase  the  resolution  for  a  given  number  of  points 
is  to  break  up  a  single  global  grid  into  several  local 
grids  that  either  overlap  each  other  or  interact 
through  a  single,  coarse,  global  grid.  The  local  grids 
could  be  generated  about  individual  components  of 
the  geometry  or  in  regions  where  the  characteristic 
length  scale  is  much  smaller  than  that  of  the  global 
flow  field,  such  u  around  the  leading  or  trailing 
edge  of  an  airfoU.  in  the  boundary  layer  or  around  a 
shock.  Since  the  individual  grids  are  only  required 
to  resolve  one  or  two  length  scales,  the  task  ot  gen¬ 
erating  them  is  greatly  simplified,  allowing  points  to 
be  clustered  where  needed  without  wasting  them 


elsewhere.  Also,  since  the  grids  are  separate,  they 
do  not  have  to  be  solved  simultaneously,  only 
requiring  information  pertaining  to  one  grid  to  be 
resident  in  memory  at  a  time.  Although,  in  gen¬ 
eral.  more  points  will  be  used  overall,  they  will  be 
used  more  efficiently  and  will  give  greater  accuracy 
for  a  given  amount  of  computing  effort. 

It  is  also  possible  to  use  different  solvers  on  the 
different  grids.  Therefore,  the  solver  that  is  most 
appropriate  to  the  purpose  of  each  grid  can  be 
chosen.  For  example,  a  full-potential  code  could  be 
solved  on  a  coarse  global  grid  while  the  Euler  equa¬ 
tions  are  solved  around  a  shock  to  correctly  predict 
the  shock  jumps  and  the  .Navier-Stokes  ^nations 
are  solved  in  the  region  near  the  body  to  resolve 
the  viscous  effects  that  are  important  there.  In  this 
way,  a  solver  that  is  only  as  sophisticated  as  it  needs 
to  be  can  be  used.  This  approach  is  similar  to  the 
singular  perturbation  techniques  in  applied 
mathematics. 

Of  course,  using  a  multiple  mesh  scheme  adds  prob¬ 
lems  inherent  to  the  approach.  Information  must 
be  passed  between  the  grids  in  such  a  way  that  the 
overall  accuracy  and  stability  of  the  scheme  is  not 
compromised  in  any  way.  Also,  the  bookkeeping 
required  to  keep  track  of  the  solution  on  several 
dinerent  grids  is  more  involved  and  needs  to  be 
automated  or  it  will  become  burdensome.  There  is 
some  memory  and  CPU  time  overhead  due  to  stor¬ 
ing  bookkeeping  information  and  setting  up  the  grid 
interaction  algorithms.  Another  possible  problem 
area  is  conservation.  Most  solvers  used  today  in 
transonic  applications  are  in  conservative  or  diver¬ 
gence  form.  This  form  is  required  for  the  scheme 
to  conserve  such  quantities  as  mass,  momentum, 
and  energy.  In  regions  of  the  flow  field  where  the 
solution  is  smooth,  conservative  form  is  not  impor¬ 
tant,  but  when  discontinuities  such  as  shocks  appear 
in  the  solution,  conservation  must  be  maintained  in 
order  to  correctly  predict  their  location.  Even  if  the 
solvers  themselves  are  conservative,  if  shocks  cross 
a  grid  boundary  that  is  not  treated  conservatively, 
the  accuracy  of  the  solution  could  suffer. 

The  method  described  here  will  interact  two  or 
more  grids  together  in  the  solution  of  unsteady, 
transonic,  viscous  flow  over  an  airfoil.  The  method 
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will  support  the  use  of  a  different  solver  for  each 
grid,  which  will  be  interacted  through  a  global 
coarse  grid  using  the  method  proposed  by  Fung  et 
aKl984).  This  approach  treats  the  local  grid  solu¬ 
tions  u  more  accurate  approaimations  to  the  correct 
flow  fleld  than  the  glob^  solution  and  uses  them  to 
correct  the  global  solution.  In  this  way.  the  global 
solution  feeu  the  effect  of  each  local  grad  and  serves 
to  transmit  this  data  to  each  of  the  other  local  grids. 
Section  2  reviews  the  approach  due  to  Fung  et  al 
and  Action  3  extends  this  approach  to  account  for 
diJCferences  in  operators  and  formulations.  In  Sec¬ 
tion  4,  we  present  results  from  three  applications  of 
the  method. 

TRUNCATION  ERROR  AND 
GRID  REFINEMENT 

Most  numerical  analysis  books  state  that  a 
differentiai  operation  L  operating  on  a  function  p  is 
related  to  a  difference  operator  (with  the  sub¬ 
script  h  denoting  the  grid  siae),  operating  on  p  by 
the  truncation  error  (abbreviated  TE  hereafter),  i.e.. 

Lp  •  L^p  ^  TE(p.k)  (1) 

This  result  of  a  Taylor  series  expansion  is  the  basis 
for  all  finite  difference  techniques.  Ordinarily,  the 
direct  system, 

L„Pk  -  0  (2) 

is  solved  for  As  a  consistency  requirement,  the 
TE  vanishes  as  the  step  size  h  approaches  zero, 
leading  to  the  limiting  solution  p 

^Um  bs  ■  d  (3) 

that  satisfies  the  differential  equation 

tb-O.  (4) 

Hence,  it  is  assumed  that  if  the  TE  is  uniformly 

small,  solving  the  discrete  system  Eqn.  (2)  will  lead 
to  a  good  approximate  solution  ba,  of  b-  However, 
in  many  p'-fl»«tems  the  TE  is  a  rapidly  varying  func¬ 
tion  of  its  arguments.  The  idea  of  convenuonal 
grid-adapting  techniques  is  to  look  for  or  modify  the 
distribution  of  g^  points  according  to  some  preset 
criteria  which  will  render  uniformity  of  TE  across 
the  solution  domain.  Unfonunately.  there  is  no 
simple  way  to  generate  a  grid  that  minimizes  the  TE 
for  a  given  problem  and,  in  many  cases,  the  process 
of  finding  the  optimal  grid  is  more  complicated  and 
time  consuming  than  computing  the  solution  itself. 

We  must  remember  that  it  is  Eqn.  (4)  that  one 
wants  to  solve,  not  Eqn.  (2).  The  discrete  equation 
that  ought  to  be  solved  corresponding  to  solving 
Eqn.  (4)  is  implied  by  Eqn.  (1),  i.e.. 

Li,p„  +  r£(b,h)  ■  0  (5) 

Here,  wejiave  deliberately  denoted  the  argument 
function  b  of  the  TE  with  a  bar,  which  can  be 
different  from  the  solution  bs-  Notice  that  if  the 
exact  solution  were  availably,  it  would  satisfy  Eqn. 
(S)  exactly,  with  b  *  bs  *  b.  This  implies  that  the 
TE  can  be  computed  exactly  by  applying  the  opera- 


tor  Li,  to  the  solution,  e.g.. 

TEfP.h)  -  -  L„p 

(fi) 

and  that  solving  the  system 

I-*b*  ■  -  TE(p,h)  ■  L„p 

(7) 

will  yield  the  exact  solution  at  nodal  points.  Hence, 
it  is  clear  that  if  the  goal  is  to  improve  a  numerical 
solution,  the  basic  grid  structure  need  not  be 
changed,  only  unproved  values  of  TE  at  grid  points 
need  to  be  provided.  To  emphasize  this  point,  the 
base  grid  in  the  examples  considered  here  is  never 
changed.  The  strategy  one  would  use  in  making 
tradeoffs  between  the  base  grid  and  approximating 
the  TE  is  not  discussed  in  this  paper. 

Analytically,  the  TE  consists  of  all  higher  deriva¬ 
tives  of  the  function  being  expanded  in  a  Taylor 
series,  if  they  exist.  If  more  neighboring  values  of 
a  fitnctiott  ue  known,  higher  derivatives  can  be 
computed,  and  hence  the  TE  can  be  better  approxi¬ 
mated.  A  TE  sequence  may  be  defined  as  follows: 

■  3’^(b*//v.h)  *  -  f'abs/v  (8) 

where  the  subscript  h/N  refers  to  values  based  on  a 
grid  of  size  h/N  (e.g.,  subdividing  the  hue  grid  of 
size  h,  N  times),  ror  simplicity,  we  may  assume 
bs//v  satisfies  the  equation 

^klN^hiN  ■  ®  (9) 

With  these  definitions  and  the  substitution  of  Eqn. 
(8)  into  Eqn.  (6)  and  setting  b  ■  bk//v>  •(  >*  easy  to 
see  that  bs  ■  bs//v  >>  a  solution  of  Eqn.  (3)  at  coin¬ 
ciding  nodal  points  (or  through  the  use  of  an  inter¬ 
polation  function). 

All  we  have  shown  so  far  is  that  it  is  possible  to 
obtain  a  refined  numerical  solution  satisfying  Eqn. 
p)  without  changing  the  base  grid,  provided  the  TE 
is  known  to  the  same  order  of  accuracy  as  the  solu¬ 
tion.  The  inclusion  of  TE  into  the  difference  equa¬ 
tion  has  been  suggested  before.  Lentini  and 
Pereyri  (1974)  proposed  a  deferred  correction  pro¬ 
cedure  to  compute  the  TE  progressively.  Warming 
and  Hyett  (1974)  and  Klopfer  and  McRae  (1983) 
implemented  it  in  forms  of  a  modified  equation. 
For  some  simple  linear  differential  equations,  even 
analytical  expressions  of  the  TE  in  terms  of  lower 
derivatives  of  the  unknown  function  have  been 
used.  However,  the  complexity  involved  in  the 
derivation  of  such  terms  and,  in  some,  the  numeri¬ 
cal  instabilities  caused  by  the  presence  of  certain 
terms  has  discouraged  the  more  popular  use  of  such 
schemes.  Liang  and  Fung  (I98S)  and  Fung  et  al 
(1984)  both  demonstrated  that  the  above  logic 
significantly  reduced  the  computation  work  required 
to  achieve  a  refined  numerical  solution  satisfying 
Eqn.  (5)  without  significantly  increasing  the  com¬ 
plexity  involved. 

TRUNCATIO.S  ERROR  DUE  TO 
OPERATOR  DIFFERENCES 

Now  extend  the  TE{P.h)  term  in  the  above  equa¬ 
tions  to  include  differences  between  the  coarse  grid 
operator  and  the  fine  grid  operator,  i.e.  the  fine 
grid  equations  can  be  a  higher  order  accuracy  ver¬ 
sion  of  the  original  equations  or  even  a  different 
equation  set  entirely.  For  example,  apply  the  above 
logic  to  a  set  of  equations  where  the  operator  can  be 
split  into  different  parts,  each  part  representing  its 
own  physics  with  its  own  characteristic  length  scales, 
e.g.  the  .Navier-Stokes  equations  (NSE).  Write  the 
NSE  in  operator  form  as 

f.,b*  +  L»b*  -  Re*’f.»''b*  *0  (10) 

where  the  total  operator  has  been  split  into  the  tem- 
poraL  convective,  and  viscous  operators.  Each  of 
these  operators  can  be  associated  with  different 
aspects  of  the  fiow  field  having  distinct  length 


IMS 


Proc.  6ih  fill-  Catrf.  on  WatJumatkai  \fodtUing 


scales.  When  these  equatioas  are  ,lK>n(  solved 
aamerically.  a  {rid  must  be  used  which  caa  capture 
aU  of  the  relevant  length  scales,  however  disparate 
they  are.  Generating  such  a  grid  which  wiu  also 
allow  the  eauations  to  be  efBciently  solved  can  be 
quite  dincnlt.  As  mentioned  earlier,  we  cu  solve 
each  operator  on  a  grid  which  resolves-  its  own 
relevant  p^ics  and  use  TE  injeaion  to  bring  the 
separate  enecu  to  a  global  grid.  For  the  sake  of 
example,  say  that  the  physical  phenomena  associ¬ 
ated  with  each  of  the  operators  m  Eqn.  ( 10)  have 
length  scales  which  decrease  from  left  to  right  and 
that  the  last  operator  is  only  important  in  some  par- 
ticulv  region  of  the  Sow  Seld.  Then,  we  can  solve 
Eqn.  (11)  on  a  grid  local  to  that  region  which  has 
sufficient  resolution  to  capture  the  relevant  physics. 

~  *0  (11) 

The  inSuence  of  this  solution  can  then  be  forced 
into  the  global  solution  via  TE  injection  by  solving 

(12) 

on  the  global  grid  which  is  only  fine  enough  to 
resolve  the  physics  relevant  to  the  L'  operator,  Le. 
M  <  M.  If  the  grids  used  in  solving  Eqn.  (11)  and 
Eqn.  (12)  cover  the  same  domain,  then  the  solu¬ 
tion  of  Eqn.  ( 12)  will  represent  the  fine  grid  solu¬ 
tion  at  the  couse  grid  nodes,  while,  if  the  fine  grui 
is  only  a  subset  of  the  total  domain,  the  solution  of 
Eqn.  (12)  will  reproduce  the  fine  grid  solution  at 
the  coarse  grid  nodes  which  lie  within  the  fine  grid 
region  with  some  influence  from  the  rest  of  the 
domain. 

Now,  if  the  temporal  effects  can  be  thought  of  as 
perturbations  to  a  steady  state,  Eqn.  (13)  can  be 
solved  for  the  final  solution  on  a  grid  which  is  only 
u  fine  u  necessary  to  resolve  the  unsteady  physics. 

(12) 

As  explained  in  Fung  et  al  ( 198T),  the  RHS  of  Eqn. 
(13)  is  an  approximation  to  the  TE  resulting  from  a 
Taylor's  senes  expansion  of  the  differential  equation 
thu  Eqn.  (12)  is  modeling  Eqn.  (13),  if  allowed  to 
converge  to  a  steady  state,  will  reproduce  the  fine 
g^  solution  at  the  coarse  grid  nodes.  For  unsteady 
calculations,  the  RHS  fixes  the  steady  solution  so 
that  only  the  unsteady  perturbations  to  this  solution 
must  be  resolved. 

In  this  section  we  have  introduced  a  procedure  simi- 
lu  to  that  suggested  by  Brandt  ( 1980)  for  separating 
the  length  scales  of  a  problem  so  that  they  can  be 
solved  more  efflciently.  Following  are  some  exam¬ 
ples  which  show  that  the  injection  of  the  TE  due  to 
differences  in  grids  and  opuators  is  a  simple  but 
effective  meant  of  improving  the  accuracy  and 
efficiency  of  a  numerical  solution. 

RESULTS 

The  first  problem  was  chosen  to  study  this  pro¬ 
cedure  with  rotated  grids  aligned  with  discontinui¬ 
ties  in  the  flow  field.  This  is  a  case  where  the  fine 
grid  operator  is  merely  a  different  form  of  the 
coarse  grid  operator.  We  solved  the  linear,  two- 
dimensional  convection-diffusion  equation  for  a 
nominal  quantity  T,  modeling  two  adjacent  fluids  of 
initially  different  temperatures  moving  at  the  same 
speed.  Upwind  differencing  was  used  for  the  con¬ 
vective  terms.  It  is  well  known  that  the  artificial 
cross-wind  difftision  due  to  upwind  differencing  is  a 
major  source  of  error.  It  causes  excessive  spreading 
of  the  discontinuity.  Fig.  1  shows  the  temperature 
contours  solved  on  a  60x40  base  grid  of  step  size  h 
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FIG.  1.  Numerical  Solution  of  a  temperature 
diffusion  layer  (60x40  grid),  no  TE  injection. 


«  0.23.  Solving  the  equations  on  a  120x80  grid 
only  slightly  improved  the  solution  and  neither  case 
compared  well  with  the  exact  solution.  Given  suit¬ 
able  pattern  recognition  schemes,  it  would  be 
natural  to  introduce  a  rotated  grid  parallel  to  the 
flow  direction  over  a  small  region  surrounding  the 
discontinuity  with  boundary  conditions  extracted 
from  the  original  solution.  Here,  we  have  done  this 
manually  with  a  20x40  grid  aligned  with  the  flow. 
Due  to  grid  rotation,  the  cross-wind  artificial 
diffusion  is  minimixed,  resulting  in  a  sharp  tempera¬ 
ture  gradient  very  close  to  the  exact  solution.  The 
isotherms  that  appear  near  the  upper  and  lower 
boundary  of  the  refined  local  solution  on  the 
subgrid  (Fig.  2)  are  an  effect  of  the  incorrect  boun¬ 
dary  conditions  extracted  from  the  base  grid  solu¬ 
tion;  these  can  be  avoided  simply  by  taking  a  larger 
subgrid.  However,  with  the  injected  TE,  the 
improved  base  solution  provides  a  sharp  gradient 
without  such  isotherms  (Fig.  3).  which  shows  that 
the  base  grid  solution  is  readjusted  smoothly 
through  the  artificial  boundaries. 


X 


FIG.  2.  Locally  refined  solution  on 
a  rotated  20x40  grid. 


The  second  example  is  the  use  of  TE  injection  to 
maintain  the  inviscid  steady  flow  corresponding  to 
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FIG.  3.  Improved  solution  of  a  temperature 
difAtsion  layer  with  TE  injection,  60x40  grid. 


body  shape  while  the  perturbed  unsteady  flow 
corresponding  to  unsteady  body  motions  are  being 
computed.  In  unsteady  transonic  flow,  for  low  to 
ffloderaM  reduced  frequencies,  the  typical  wave 
length  of  an  acoustic  signal  is  of  the  order  of  the 
chord  of  the  wing.  Hence,  a  grid  with  a  minimum 
spacing  of  a  tenth  of  the  chord  should  be 
suf&iently  flne  to  resolve  these  waves.  However, 
In  a  typical  inviscid  flow,  the  smallest  characteristic 
length  scale  is  on  the  order  of  the  radius  of  curva¬ 
ture  of  Che  leading  edge  of  the  airfoiL  Therefore,  a 
grid  with  a  minimum  spacing  of  a  hundredth  of  a 
chord  is  needed  to  resolve  the  flow  Held.  Due  to 
linear  or  nonlinear  numerical  instability,  these  grids 
may  require  too  small  a  time  step  for  efflcient  com¬ 
putations  of  the  unsteady  acoustic  waves  caused  by 
the  small  unsteady  wing  motioos  and  deformation 
assumed  in  flutter  analysis,  la  order  to  ease  the 
restrictions  on  the  time  step,  Fung  and  Fu(198S) 
used  the  technique  described  above  to  compute  the 
steady  and  unsteady  flows  on  different  grids. 

They  usumed  inviscid  flow  so  Re-w  and  only 
Eqns.  ( 12)  and  ( 13)  above  are  relevant.  Hqn.  (12) 
is  solved  on  a  pH  which  is  line  enough  to  resolve 
all  the  relevant  physics  in  the  steady  flow  field. 
Eqn.  (13)  is  solved  on  a  grid  which  is  only  as  flne  as 
required  to  resolve  the  acoustic  signals  at  dcKribed 
above.  Results  are  presented  for  supercritical  flow 
over  an  NLR  7301  airfoil.  The  airfoil  is  pitching 
harmonically  at  a  reduced  frequency  k  ■  6. 1  in  a 
freestream  with  Mach  number  of  0.7.  The  calcula¬ 
tions  were  performed  on  a  fine  grid  109x87,  a 
coarse  grid  54x43,  and  on  the  coarse  grid  with  TE 
injection  from  the  fine  gri^  Fig.  4  compares  the 
unsteady  pressure  distributions  at  90  degree  inter¬ 
vals.  From  these  comparisons,  it  it  evident  that, 
except  for  minor  differences  near  the  shock,  the 
resnlu  obtained  on  the  coarse  grid  with  TE  injection 
are  just  u  accurate  as  those  on  the  fine  grid.  The 
fine  grid  solution  required  900  seconds  of  CPU 
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FIG.  4.  Comparison  of  the  unsteady  pressure 
distributiont  of  an  NLR  7301  atrfoU. 


time,  the  coarse  grid  solution  57  seconds,  and  the 
coarse  grid  solution  with  TE  injection  60  seconds. 
A  reduction  in  computation  time  by  a  factor  of  four 
is  attributed  to  the  reduction  in  grid  points  and 
another  factor  of  four  to  the  relaxation  of  the  allow¬ 
able  time  step  imposed  by  numerical  stability. 

This  technique  has  also  been  applied  by  Schoen 
(1987)  to  three  dimensional  flow  in  a  code 
developed  at  Nasa-Ames  for  Transonic  UNsteady 
Aerodynamics  •  TUNA,  due  to  Bridgeman  and 
Sieger  (1982).  Results  are  presented  for  supercnti- 
cai  flow  over  a  rectangular  wing  with  NACA-(X)12 
cross-section  and  an  aspect  ratio  of  6.  Fig.  5  shows 
the  deviation  of  the  pressure  coefficient  from  the 
fine  grid  stesdy  state  at  the  airfoil  midspan  for  the 
fine  gnd  solution  and  the  coarse  grid  solution  with 
TE  injection.  Note  that  the  scales  used  in  Fig.  5  do 
not  permit  us  to  show  the  large  deviations  found  in 
the  coarse  grid  solution  without  TE  injection.  The 
wing  is  plunging  harmonically  at  a  reduced  fre¬ 
quency  of  0.4  in  a  freestream  with  Mach  number 
0.7.  The  calculations  were  performed  on  a  fine  grid 
89x49x18,  on  a  coarse  grid  45x25x18,  and  on  the 
coarse  grid  with  TE  injected.  The  grids  were  not 
coarsened  in  the  spanwise  direction  as  the  typical 
spacing  in  this  dimension  is  too  large  on  even  the 
fine  grid  for  resolving  acoustic  waves.  The  method 
captures  the  steady-state  solution  exactly.  Various 
phases  from  the  first  cycle  are  shown  to  verify  that 
the  unsteady  procedure  maintains  respectable  accu¬ 
racy  u  well,  with  100  steps  per  cycle,  the  fine  grid 
solution  took  68  Cray  seconds  per  cycle  and  the 
coarse  grid  with  or  without  TE  took  19  seconds  per 
cycle.  Because  the  time  marching  scheme  in  this 
code  is  unconditionally  stable,  there  were  no  time 
step  restrictions  to  be  relaxed  on  the  coarse  grid. 
Thus,  a  reduction  in  compulation  time  by  a  factor 
of  about  four  is  attributed  to  the  reduction  in  grid 
points. 

The  third  example  shows  the  use  of  TE  injection  to 
fix  the  viscous  effects  in  a  steady  flow  field  while 
using  a  less  dense  grid  to  calculate  the  inviscid 
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FIG.  5.  Comparisons  of  the  unsteady 
deviations  from  the  fine  grid  steady  state 
solution  for  an  ONERA-M6  airfoU. 


aspects  of  the  flow.  As  mentioned  above,  a  grid 
with  a  minimum  spacing  of  a  hundredth  of  a  chord 
should  be  sufficient  to  resolve  the  inviscid  flow 
field.  However,  in  a  viscous  calculation,  the  impor* 
tant  length  scale  is  thgi^f  the  boundary  layer  thick' 
ness  which  is  0(l/vRe),  and  in  a  turbulent  flow 
there  is  a  sublayer  whose  length  scale  is  even 
smaller.  For  accuracy,  a  grid  must  be  used  which 
will  resolve  these  smalt  scales  and  yet  must  still 
retain  enough  points  away  from  the  body  to  resolve 
the  far  field.  These  requirements,  coupled  with 
smoothness  limitations,  result  in  a  highly  clustered 
grid  with  a  large  number  of  nodes  which  tends  to 
slow  convergence  of  viscous  calculations.  In  order 
to  resolve  the  two  regions  of  flow  with  their  distinct 
length  scales  more  efficiently,  we  have  applied  the 
above  method  to  allow  she  viscous  and  inviscid  cal¬ 
culations  to  be  computed  on  separate  grids. 

The  method  is  implemented  using  the  thin-layer 
N avier- Stokes/ Eli ler  code  ARC20  due  to  Steger 
(1978)  and  extended  by  Pulliam  and  Chaussee 
(1981,1984).  The  code  was  applied  to  a  laminar, 
compressible  boundary  layer  on  a  flat  plate  at  zero 
angle  of  attack;  see  White  ( 1974)  for  details  of  this 
solution. 

A  laminar,  viscous  solution  was  found  on  a  61x61 
base  grid,  corresponding  to  the  solution  of  Eqn. 
(11),  m  the  full  domain  using  ARC20  and  desig¬ 
nated  the  fine  grid  solution.  At  first,  in  Eqn.  (12). 
the  forcing  terms  were  formed  using  the  tame  grid 
as  the  coarse  grid.  The  forced  calculation  was  then 
performed  with  free  stream  initial  conditions.  As 
expected,  the  forced  calculation  returned  the  origi¬ 
nal  solution. 

Eight  more  grids  were  formed  from  the  base  grid  by 
keeping  every  other  or  every  fourth  coordinate  line 
in  either  direction.  The  fine  grid  solution  was  res¬ 
tricted  to  each  of  the  eight  coarser  grids.  These 
eight  grids  are  subsets  of  the  fine  grid  so  the  restric¬ 
tion  process  is  merely  an  injection  of  the  fine  grid 
solution  at  common  grid  points.  In  each  case,  the 
forcing  function  was  formed  and  a  forced  calcula¬ 
tion  performed  using  free  stream  initial  conditions. 
Only  the  results  from  the  finest  and  coarsest  grids 
will  be  shown  here. 
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FIG.  6.  Fine  (61x61)  and  coarse  (16x16)  grids 
used  for  calculation  of  boundary  layer  flow. 


Fig.  6  shows  the  two  grids  used.  Fig.  7  compares 
the  results  found  on  the  16x16  grid  using  the  Euler 
equations  plus  forcing  terms  (INVISJ)  with  the  fine 
grid  results  using  thin-layer  Navier-Siokes  (V1S4). 
As  the  plot  shows,  the  horizontal  momentum  is 
resolved  very  well  by  the  coarse  grid  when  forcing 
terms  are  used.  Fig.  7  also  compares  INVISJ  with 
thin-layer  Navier-Stokes  results  on  the  16x16  grid 
without  forcing(VlSD).  It  is  obvious  that  INVISJ  is 
a  far  better  solution  than  VISD. 

CONCLUDING  RE.MARKS 

The  technique  of  truncation  error  injection  intro¬ 
duced  by  Fung  et  al  (1987)  has  been  successfully 
generalized  to  include  operator  differences  in  the 
truncation  error  term.  This  generalization  has 
enabled  the  use  of  multiple  grids  to  resolve  widely 
disparate  length  scales,  signmcantly  increasing  the 
efficiency  of  the  solution  process  while  retaining 
accuracy.  The  different  grids  are  each  required  to 
resolve  only  one  aspect  of  the  flow  field,  making 
them  much  easier  to  generate  and  more  efficient  in 
their  use  of  points. 
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FIG.  7.  Comparison  of  coarse  grid  Euler  with 
forcing  results  (INVISJ)  with  fine  grid 
Navier-Stokes  results  (VIS4),  ^nd  coarse  grid 
Navier  Stokes  results  (VISD). 
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IntrodiKtion 

HE  urgent  need  for  effective,  reliable  methods  for 
unsteady  aerodynamic  predictions  at  transonic  Mach 
numbers  is  evident  from  the  Fanner  and  Hanson'  experiment 
in  which  it  was  observed  that  the  flutter  boundary  for  a  wing 
with  a  supercritical  cross  section  is  substantially  lower  than 
that  with  a  conventional  one.  At  transonic  speeds,  the  size  and 
location  of  the  embedded  supersonic  zone  over  the  wing  affect 
the  way  acoustic  signals  propagate  and.  hence,  the  aerody¬ 
namic  responses  to  disturbances.  Recent  developments  in 
computational  fluid  dynamics  and  the  availability  of  super¬ 
computers  have  made  accurate  flow  prediction  possible. 
However,  for  applications  like  routine  flutter  calculation  and 
aircraft  design  optimization,  the  currently  available  codes, 
especially  the  ones  for  three-dimensional  computations  like 
XTRAN3S  of  Rizzetta  and  Borland’  and  USTF3  of  Isogai  and 
Suetsugu,’  are  still  much  too  time  consuming. 

As  mentioned  in  Fung.*  one  of  the  problems  in  unsteady 
transonic  flow  computation  is  the  grid  for  obtaining  the  solu¬ 
tion.  Aside  from  the  issue  of  finding  the  best  grid  for  a  given 
wing  geometry,  a  grid  mutt  have  a  local  mesh  size  comparable 
to  the  radius  of  curvature  of  the  leading  edge  in  order  to  prop¬ 
erly  resolve  the  fast  expansion  that  determines  the  size  of  the 
sonic  bubble  and  the  strength  and  location  of  the  shock.  The 
computational  domain  must  be  large  enough  to  allow  the  flow 
to  relax  to  the  freestream  condition  with  little  confinement 
from  grid  boundaries.  Almost  all  grids  currently  used  for 
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aerodynamic  compulations  are  based  on  these  considerations. 
However,  these  grids,  while  suitable  for  computing  steady 
flows,  may  require  (due  to  linear  or  nonlinear  numerical  in¬ 
stability)  too  small  a  time-step  limitation  for  efficient  com¬ 
pulations  of  the  unsteady  acoustic  waves  due  to  the  small 
unsteady  wing  motions  and  deformations  assumed  in  flutter 
analysis.  For  low-to-moderate  reduced  frequencies,  the  typical 
wavelength  of  an  acoustic  signal  in  a  transonic  flow  is  of  the 
order  of  the  chord  of  the  wing.  Hence,  a  grid  with  a  minimum 
spacing  of  a  tenth  of  the  chord  should  be  sufficient.  However, 
an  accurate  prediction  of  the  steady  flowfield  over  a  wing  at 
supercritical  Mach  number  often  requires  a  minimum  spacing 
of  a  hundredth  of  the  chord  and.  hence,  a  time-step  require¬ 
ment  based  on  the  CFL  condition  10  times  as  restricted  as  that 
needed  for  accuracy. 

In  this  Note,  a  technique  is  introduced  that  allows  the  steady 
and  unsteady  flows  to  be  computed  on  different  grids.  To 
demonstrate  the  efficiency  of  this  technique,  the  unsteady 
small-disturbance  transonic  equation  is  used  for  unsteady 
aerodynamic  prediction.  The  results  of  applying  this  technique 
are  compared  to  those  obtained  on  single  grids. 

Computaiions  with  TmneatioR  Error  Injection 

It  has  been  shown’  that,  by  solving  the  corresponding  dif¬ 
ference  equation  with  the  truncation  error  includinl  as  a  forc¬ 
ing  term,  exact  nodal  values  of  the  solution  to  a  differential 
equation  can  be  obtained;  that  if  the  exact  solution  were 
known,  the  exact  truncation  error  can  be  computed  at  nodal 
points;  and  that  the  truncation  error  can  be  approximated  by 
local  grid  refinement. 

Consider  a  difference  equation  of  the  form 

(1) 

where  L,  and  L«  correspond  to  the  temporal  and  spatial 
discrete  operators,  respectively,  and  the  numerical  solution 
on  a  grid  of  size  h.  Assuming  that  is  the  steady-state  solu¬ 
tion  of  Eq.  (I)  satisfying 

LM-0  (2) 

it  is  quite  obvious  that  also  satisfies 

(3) 

and  that  solving  Eq.  (3)  for  0,  is  the  same  as  solving  Eq.  (I). 
However.  Eq.  (3)  is  more  general  in  the  sense  that,  if  it  con¬ 
verges,  it  will  yield  the  steady  state  oi  regardless  of  whether  oj 
satisfies  Eq.  (2).  For  example,  we  could  replace  oj  by 
i.e., 

*  f-»d»  ■  f.»dil/;v  (^) 
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Fig.  I  Cemparttaa  of  ibe  tarfaet  prtsiurt  disirtbaitow  of  aa  NLR 
T3tl  alrfotl  M  <•  O.TV  oMaiatd  ^  dirrtfcM  awtbodi  oa  dirfcteai 
grids. 
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TiM*  I  Comptrbom  af  cawpatMiaa  liaw  far  dirfmai 
_ caapamlaai.  CPD’y _ 


NACA-»»A-<?1Q  NCR  7}0l 


Case 

Steady 

Unsteady 

Total 

Steady 

Unsteady 

Total 

Fine  grid 

13.2 

299.1 

312.3 

32.9 

«99  5 

932.4 

Coarse  grid 

3.J 

14.3 

20.0 

10.4 

56.3 

66.9 

Coarse  ♦  TEI 

17.0 

l«.) 

33.3 

36.7 

59.3 

96.0 

Here,  oS/.v  denotes  the  steady-state  solution  on  a  finer  grid  of 
size  b/yv,  satisfying  the  difference  equation 

(J) 

As  explained  in  Fung  et  al.,’  the  term  in  Eq.  (4)  is  an 

approximation  to  the  truncation  error  resulting  from  a 
Taylor’s  series  expansion  of  the  differential  equation  that  Eq. 
(1)  is  modeling,  liie  difference  between  Eqs.  (3)  and  (4)  is  that 
(he  latter  contains  an  approximation  of  the  truncation  error  of 
the  steady  solution  due  to  discretization  and  hence  will  yield 
the  more  accurate  as  a  steady  solution  if,  of  course,  the 
unsteady  effects  are  allowed  to  diminish.  Assuming  that  the 
base  grid  of  size  A  is  fine  enough  for  resolving  the  unsteady 
part  of  there  will  be  very  little  or  no  truncation  error  due 
to  discretization  and  Eq.  (4)  should  yield  more  accurate  solu¬ 
tions  that  Eq.  (I). 

In  the  case  where  the  spatial  operator  is  linear,  Eq.  (4)  is 
simply  Eq  (I)  with  the  term  L„^X,s  added  to  both  sides  of  it. 
For  nonlinear  operators,  Eq.  (4)  will  enforce  the  penurbations 
about  a  steady  state  to  satisfy  the  governing  equation.  It  is  of 
interest  to  note  that  the  splitting  of  an  operator  into  a  steady 
and  an  unsteady  operator  as  in  Eq.  (1)  is  necessary  only  for  the 
clarity  of  the  discussion. 
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ng.  3  Real  aaO  laiaglaary  parts  of  Ibe  aasMady  pntsarc  bisiribaiioa 
at  aa  NLR  73*1  lUrtell  pMcMag  baraMUtcally  at  M.  -II.TO  aad  at 
reduesd  trtqaaacy  b ••.!  obtaiaaO  by  diffetaat  aMboea  oa  tUtrereai 
grtda. 

NuincfteRl  Example 

To  demonstrate  its  effectiveness,  we  have  implemented  this 
technique  into  the  code  AZTRAN*  for  solving  the  unsteady 
transonic  small-disturbance  equation.  No  modification  to  the 
code  is  needed  except  adding  the  stored  truncation  error  values 
computed  from  the  steady-state  solution  to  the  matrix  equa¬ 
tions  at  each  time  step.  Figure  I  shows  different  steady-state 
pressure  distributions  for  an  NLR  7301  airfoil  at  a  freestream 
Mach  number  Af.  -0.70  obtained  on  a  fine  grid  (109  x  87) 
and  on  a  coarse  grid  (34  x  43)  that  was  formed  by  omitting 
every  other  grid  line  of  the  fine  grid.  As  a  result  of  faster  ex¬ 
pansion  at  the  leading  edge,  the  pressure  distribution  (solid 
line)  obtained  on  the  finer  grid  shows  a  stronser  shock  located 
further  downstream  compared  to  that  (triangles)  obtained  on 
the  coarser  grid.  The  pressure  distribution  shown  by  the 
squares  was  obtained  on  the  coarse  grid  with  injected  trunca¬ 
tion  error  computed  from  the  fine-grid  solution.  It  coincides 
with  the  solid  line  almost  everywhere  except  near  the  shock, 
where  a  smearing  of  the  preuure  jump  occurs  mainly  because 
of  numerical  differentiation  and  graphic  interpolation.  The 
numerical  smearing  of  the  shock  within  one  grid  point, 
however,  is  of  secondary  importance  as  far  as  the  computation 
of  integrated  loads  for  an  aeroelastk  application  is  concerned. 
Theoretically,  unsteady  penurbations  on  a  shock  create  a 
singularity,  an  integrable  one,  however,  in  the  penurbed 
quantities  at  the  mean  shock  position.  In  reality,  this  singular¬ 
ity  is  often  smeared  by  viscous  effects.  Motions  of  the  shock,  a 
nonlinear  phenomenon  in  nature,  also  cause  higher  harmonics 
in  the  aerodynamic  forces.  However,  it  was  shown  by  Davis 
and  Malcolm*  that  the  higher  harmonics  have  little  effect  on 
the  integrated  loads. 

Corresponding  unsteady  pressure  distributions  at  90-deg  in¬ 
tervals  for  the  same  airfoil  pitching  harmonically  at  a  reduced 
frequency  A>i0.l  are  shown  in  Fig.  2  and  their  Hrst  harmonic 
decompositions  in  Fig.  3.  It  is  evident  from  these  comparisons 
that,  except  for  the  minor  differences  near  the  shock,  the 
results,  including  the  motions  of  the  shock,  obtained  on  the 
coarse  grid  with  truncation  error  injection  are  just  as  accurate 
as  those  on  the  fine  grid.  The  CPU  time,  as  listed  in  Table  I, 
for  obtaining  five  cycles  of  harmonic  motion  on  the  fine  grid 
was  900  s,  for  the  coarse  grid  solution  it  was  37  s,  and  for  the 
coarse  grid  solution  with  the  truncation  error  injection  it  was 
60  s.  A  reduction  in  computation  time  by  a  factor  of  four  is  at¬ 
tributed  to  the  reduction  in  grid  points  and  another  factor  of 
four  to  the  relaxation  of  the  allowable  time  step  imposed  by 
numerical  stability. 

Similar  results  for  a  conventional  airfoil,  NACA-64A-0I0, 
at  a  freestream  Mach  number  of  0.796  and  a  reduced  fre¬ 
quency  of  0. 1  are  also  listed  in  Table  I .  Because  of  the  conven¬ 
tional  leading-edge  shape,  the  pressure  expansion  at  the 
leading  edge  is  quite  miM  compared  to  that  of  the  NLR  7301 . 
Hence,  the  results  obtained  on  both  grids  were  just  as  ac¬ 
curate,  except  near  the  shock,  with  or  without  the  truncation 
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error  injection.  This  verifies  an  assumption  in  our  theory  that 
a  rciaciveiy  coarse  grid  is  sufficient  for  resolving  the  unsteady 
acoustic  waves  if  the  steady  pressure  is  accurate. 

Conclusioa 

A  simple  numerical  technique  which  can  be  easily  im¬ 
plemented  in  any  numerical  code  for  computations  of  two-  or 
three-dimensional  unsteady  transonic  aerodynamics  has  been 
introduced.  This  technique  allows  the  decoupling  of  a  solution 
having  two  distinct  length  scales  into  two  parts.  By  solving 
each  part  on  a  grid  of  the  proper  length  Kale,  substantial  sav¬ 
ings  in  computation  time  and  storage  can  be  achieved. 
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latrodoctioR 

HE  problem  of  intereu  is  flow  ac  subsonic  freestream 
Mach  number  past  a  wedge  (see  Fig.  1).  After  a 
compression  at  the  leading  edge,  the  flow  expands  and 
reaches  the  sonic  condition  at  the  shoulder.  Downstream  of 
the  shoulder,  the  flow  must  return  to  the  subsonic  freestream 
condition  through  a  shock,  or  shocks.  Because  of  its 
simplicity,  this  model  was  us^  by  Cole'  and  Yoshihara*  to 
study  the  characteristics  of  transonic  flow  and  the  validity  of 
the  transonic  small-disturbance  equation.  However,  because 
of  the  inherent  limitations  of  the  hodograph  method  used  in 
their  studies,  certain  assumptions  regarding  the  flow  pattern 
were  made  in  order  to  have  a  complete  specincation  of  the 
problem  in  the  hodograph  plane.  Cole  assumed  that  the 
sonic  line  is  locally  vertietd  and  contended  that  the  PrandtI- 
Meyer  expansion  at  the  wedge  shoulder  be  terminated  by  an 
oblique  shock,  whereas  Yoshihara  assumed  that  a  smooth 
overexpansion  occurs  at  the  shoulder  and  the  supersonic 
zone  is  terminated  by  a  normal  shock  downstream  of  the 
shoulder.  The  flow  was  studied  experimentally  by  Liepniann 
and  Bryson,’  whose  results  were  inconclusive  about  the 
oblique  shock  at  the  shoulder  due  to  viscous  effects,  and 
numerically  by  Yu  and  Seebass,*  whose  solution  was  inac¬ 
curate  due  to  the  numerical  scheme  used  and  insufficient 
resolution  at  the  shoulder.  To  date,  no  accurate  solution 
describing  the  correa  flow  structure  has  been  reported. 

Realizing  that  the  flow  structure  is  a  result  of  the  expan¬ 
sion  at  the  shoulder  where  small  variations  in  flow  variables 
may  affect  the  overall  flow  structure,  we  apply  a  subgrid 
refinement  procedure  (suggested  by  Fung  et  al.’),  which 
allows  local,  as  well  as  global,  improvement  to  obtain  a 
refined  solution  on  a  fixed  grid  without  changing  the  base 
grid  structure,  as  is  needed  in  other  grid  refinement  pro¬ 
cedures.  In  this  procedure,  a  solution  is  defined  on  a  fixed 
grid  and  is  improved  through  approximation  of  the  trunca¬ 
tion  error,  which  is  usually  ignored  in  a  conventional  finiie- 
difference  approximation.  It  is  demonstrated  in  this  study 
that  local  refinement  yields  results  as  accurate  as  those  ob¬ 
tained  on  a  uniformly  refined  grid,  that  the  truncation  error 
is  an  effective  means  for  measuring  error  in  the  solution  of  a 
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Kung  University,  Taiwan. 
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difference  equation  and  for  improving  a  numerical  solution 
on  a  fixed  grid,  and  that  the  procedure  suggested  by  Fung  et 
al.’  achieves  substantial  savings  in  computer  time  and 
storage  with  minimal  bookkeeping  effort. 

A  sequence  of  improved  solutions  obtained  using  this 
method  indicates  the  existence  of  both  an  oblique  shock  near 
the  shoulder  and  a  normal  shock  downstream  of  it. 

Governing  Equnlion  and  Numerical  Scheme 

The  governing  equation  for  this  flow  is  the  transonic, 
small-disturbance  potential  equation; 

(I) 

Here,  the  perturbed  velocity  potential  d  and  the  space  coor¬ 
dinates  X  and  y  have  been  scaled  properly.  A',  the  transonic 
similarity  parameter,  is  a  result  of  the  scaling,  which  involves 
the  freestream  Mach  number,  the  wedge  half  thickness  9, 
and  the  ratio  of  specific  heats  -y.  The  boundary  condition  on 
the  body  is  the  tangency  condition. 

d,(jr,0)*0,  jr>  I 

>1,  Osxsl 

and  the  far-field  boundary  conditions  are 

d,  »  d,  »  0  as  x’  +  y’  —  a» 

For  computational  efficiency,  this  condition  may  be  replaced 
by  an  analytical  expression  for  d  corresponding  to  a  source 
singularity  at  point  ('/>,  0). 

Equation  (I)  admits  weak  solutions  that  satisfy  the  shock 
jump  condition 

(/f-d,)[dj’  + [d,)’=0  (2) 

where  the  tilde  denotes  an  average  quantity  across  the  shock 
and  the  brackets  denote  the  jump  in  the  argument. 

A  balance  of  x  momentum  requires  that  the  pressure  drag 
on  the  wedge  be  related  to  the  shock  jumps,  i.e., 

f  C,d,(x,0)dx- -  ['jd.dx-— [O.l’dy  (3) 

J bo^Y  Jo  6  J <bocli« 

This  equation  is  used  later  for  evaluating  the  accuracy  of  the 
results. 

Equation  (I)  is  discretized  using  the  monotone  scheme  sug¬ 
gested  by  Engquist  and  Osher,‘  and  the  resultant  difference 
equations  are  solved  by  a  tine  relaxation  method. 
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Tnincalion  Error  and  Grid  RcnnemenI 

The  dirference  between  a  differential  equation  £O>0  and 
the  difference  equation  that  models  it  is  the  truncation  error 
(denoted  as  TE  hereafter),  i.e., 

£0-L*^-t-TE(0,/i) 

Here,  £  is  the  differential  operator;  p,  the  solution;  and  L«, 
the  difference  operator  resulting  from  retention  of  ap¬ 
propriate  leading  terms  in  a  Taylor’s  series  expansion  on  a 
grid  of  typical  spacing  h. 

Ordinarily,  jhe  discrete  system,  is  solved  for  an 

approximate  i,,  and  the  TE  is  assumed  to  be  small.  Unfor¬ 
tunately,  unless  the  local  grid  size  is  comparable  to,  or 
smaller  than,  the  local  length  scale  of  the  Mlution,  the  local 
TE  will  be  large,  causing  an  error  in  Attempts  and 
limited  successes  have  been  reported  on  various  adaptive  grid 
distribution  procedures  (e.g.,  Pearson’  and  Babuska*), 
which  tend  to  minimize  the  local  TE.  However,  these 
methods  very  often  depend  on  the  criteria  used  and  are  sen¬ 
sitive  to  the  control  parameters  introduced  in  the  procedure 
for  relating  the  grid  sizes  to  the  measures  of  error.  In  Ref.  5, 
Fung  et  al.  explained  that  it  should  be  possible  to  obtain  at 
the  node  points  of  a  grid  the  values  of  the  exact  solution  to  a 
differential  equation  if  the  TE  were  known.  This  means  that 
if  the  goal  is  to  improve  the  numerical  solution,  the  basic 
grid  structure  does  not  have  to  be  changed— only  improved 
values  of  the  TE  need  to  be  provided. 

It  can  be  shown  that  if  the  exact  solution  •  were  known, 
the  TE  can  be  computed  exactly  by  applying  the  operator  L, 
to  e.g., 

TE(d,A)  ■  -L*d 

Hence,  if  o,  ^  is  the  refined  solution  on  the  grid  of  size  A/V 
(i.e.,  subdividing  the  base  grid  of  size  h.  N  times),  the 
operation 

v*TE»  v»TE(^./i) 

would  yield  better  approximations  of  the  TE  as  the  refine¬ 
ment  factor  N  increases.  Of  course,  this  is  true  only  if  the 
scheme  is  consistent  with  the  differential  equation  being  ap¬ 
proximated;  i.e.,  the  truncation  error  on  the  finest  subgrid 
diminishes  while  the  truncation  error  on  the  base  grid  ap¬ 
proaches  the  asymptote. 

To  effect  an  improved  solution  on  a  Fixed  grid  by  TE  in¬ 
clusion.  we  use  the  procedure  suggested  by  Fung  et  al.’  To 
begin,  we  choose  a  base  grid,  which  will  not  change  through¬ 
out  the  procedure.  The  TE  is  set  to  zero,  and  the  discrete 
equation 

L***-i-TE-0  (4) 

is  solved  for  the  fint  approximation  Regarding  d/|  as  a 
refined  solution  for  a  grid  of  size  2A,  an  estimate  of  the  TE 
is  then  formed  by  computing 

The  regions  where  the  estimated  TE,  -  is  larger  than 
a  preset  value  r  are  identified;  if  the  base  grid  was  properly 


chosen,  these  regions  should  be  a  small  part  of  the  base 
region.  With  boundary  conditions  interpolated  from  the  base 
solution,  refined  solutions  satisfying 

L*/A(d*/AI  =  0 

are  then  obained  and  used  to  form  the  approximated  TE, 
Outside  the  region  where  refinement  is  needed, 
the  TE  is  assumed  to  remain  zero.  The  process  is  repeated 
with  the  newly  obtained  TE»,v  until 

<5  for  all 

This  procedure  is  shown  schematically  in  Fig.  2. 

It  was  shown  by  Fung  et  al.’  how,  for  various  cases,  this 
procedure  can  lead  to  improvement  of  numerical  solution  on 
a  fixed  grid,  including  the  sharpness  of  shocks.  A  similar 
procedure  has  been  suggested  by  Brandt.*  He  pointed  out 
that  nesting  subgrids  can  be  coupled  with  the  multigrid  pro¬ 
cedure,  which  has  been  widely  used  for  accelerating  the  con¬ 
vergence  of  numerical  schemes. 

Results  and  Discussion 

The  base  grid  for  all  of  the  numerical  results  is  a  uniform 
grid  of  spacing  haO.I  over  a  computational  domain  of  size 
6x6  wedge  chords.  Figure  3  shows  the  distribution  of  the 
estimated  TE,  computed  after  obtaining  the  initial 

base  solution  pj.  The  values  shown  have  been  multiplied  by 
a  factor  of  10.  Subregi>  }  outside  of  which  the  maximum 
TE  is  less  than  the  preset  t  are  shown  in  Fig.  4.  No  attempts 
were  made  to  tailor  these  regions  to  the  minimal  sizes.  The 
case  of  rwO  means  that  the  whole  domain  is  refined.  For  the 
case  of  rw3.0,  two  regions  were  identified;  one  is  about  the 
leading  edge  and  the  other  about  the  shoulder,  as  expected. 
Local  solutions  over  these  regions  were  then  obtained  with 
interpolated  boundary  conditions  from  pj  (for  details,  see 
Ref.  10).  Prom  the  local  refined  solution,  the  approximated 
truncation  error  is  then  computed  and  used  to  obtain  the 
next  refined  solution.  Pi,  on  the  base  grid  satisfying  Eq.  (4). 


YES 


ng.  I  flew  over  a  wedge  «  a  traasoak  Mach  aaiaher. 


Fig.  2  Scheanlic  of  adaptive  grM  rcflaencnt  with  iraacalloa  error 
injection. 
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TaMc  I  Cowparisoa  of  drag  for  olti  of  f  after  the  flol  reHawBcal  cycte 


Rcfiacincm 
cycle,  k 

Preset 
loleraoce.  r 

Pressure 

drag 

Wave  drae  due  to  shocks 
Oblique  Normal  Total 

CPU  on 
CRAY  X-MP.  5 

Rdative 
error  of 
wave  drag, 

0 

0.750 

0.300 

0.538 

0.838 

33 

12 

1 

0.00 

0.772 

0.201 

0.610 

0.811 

629 

5 

1 

0.05 

0.773 

0.202 

0.608 

0.810 

66 

5 

1 

0.50 

0.770 

0.179 

0.627 

0.806 

S3 

5 

1 

3.00 

0.784 

0.244 

0.575 

0.819 

50 

5 

Table  2  Coapatiaaa  of  the  drag  values  after  dirreftat  estdes  of  renucaMot, 

,  r«0.05 

Refinement  Grid  refinement 
cycle,  k  factor,  .V 

Pressure 

drag 

Wave  drat  due  to  shocks 
Oblique  Normal  Total 

CPU  on 
CRAY  X-.VIP. 

Relative 
error  of 

$  wave  drag,  rii 

0  1 

0.750 

0.300 

0.538 

0.838 

33 

12 

1  2 

0.773 

0.202 

0.608 

0.810 

66 

5 

2  2 

0.773 

0.202 

0.608 

0.810 

78 

5 

2  4 

0.793 

0.129 

0.672 

0.801 

255 

1 
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Fig.  4  Compaitaliooal  doataio  aad  labrtgioaa  for  vtrtoaa  grid 
rtfiaaowalt. 


Figure  5  compares  the  base  solutions  before  and  after  the 
nrsi  cycle  of  refinement  with  that  obtained  on  a  uniformly 
refined  grid  size  of  /i/2.  It  shows  that  this  refinement  pro¬ 
cedure  is  accurate  and  efficient.  The  effects  of  the  refine¬ 
ment  with  TE  injection  on  the  sonic  line  and  on  the  shock 
are  shown  in  Fig.  6.  Because  of  the  improvement  at  the 
shoulder,  the  normal  shock  moves  downstream  and  becomes 
stronger.  Table  I  compares  values  of  the  pressure  drag  on 
the  shoulder  and  the  wave  drag  due  to  the  shocks  (according 


to  Eq.  (3)1  obtained  for  different  values  of  r  which  corres¬ 
pond  to  different  regions  of  refinement.  It  is  quite  clear 
from  this  comparison  that  a  subregion  of  r  »  O.OS  is  adequate 
for  accurate  calculations.  The  difference  between  the  wave 
and  the  pressure  drag  diminishes  after  one  cycle  of  refine¬ 
ment  from  12  to  at  a  cost  of  only  twice  the  computation 
time  needed  for  the  base  solution.  This  is  a  substantial  sav¬ 
ings  relative  to  a  uniformly  refined  solution,  which  costs  20 
times  as  much  as  the  base  grid  solution  oj.  It  is  shown  that 
further  savings  in  computation  time  are  possible  with  minor 
losses  in  accuracy.  Another  cycle  of  refinement,  keeping  the 
refinement  factor  the  same  (At>2)  but  using  boundary  con¬ 
ditions  from  the  base  solution  for  the  local  refined  solution, 
showed  that  the  effects  due  to  the  newly  updated  boundary 
conditions  on  the  base  solution  were  minor  compared  to 
those  of  increasing  JV.  The  difference  in  the  drag  values  is 
reduced  to  less  than  I*’*  after  two  cycles  of  refinement, 
Ar«2,  and  for  a  refinement  factor  /V«4,  but  to  only  SW*  if 
A(»2  (Table  2).  For  two  cycles  with  N=4,  the  computation 
time  was  2SS  s,  approximately  four  times  that  for  one  cycle 
because  the  total  number  of  grid  points  was  increased  by  the 
same  factor:  a  better  strategy  would  be  to  apply  the  adaptive 
refinement  procedure  to  the  local  solution  coupled  with  a 
pattern  recognition  routine  to  determine  the  subgrid 
geometry,  as  demonstrated  by  Berger  and  Oliger,"  instead 
of  the  uniform  refinement.  The  base  solution  obtained  after 
two  cycles  of  refinement  was  accurate  within  1  %  in  the  drag 
values,  the  same  order  of  magnitude  as  the  maximum 
residual  of  the  difference  equation.  Further  refinements  are 
not  needed  unless  the  residual  tolerance  is  set  to  a  smaller 
value  in  solving  the  difference  equations. 
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Fit.  7  {Mstribatio*  of  momcMaa  Umsm  laiMt  a  eonlroi  toIubm  tl 
dlffcrtM  ralaat  of  y.  X^O.S. 


Fit.  •  Soak  Hat  tad  obllaat  shock,  KmQ.S. 


In  evaluating  the  shock  drag,  one  notices  that  the  drag  due 
to  the  normal  shock  is  about  four-fifths  of  the  pressure  drag 
on  the  shoulder.  To  detect  the  momentum  loss,  we  moved  a 
control  volume  along  a  grid  line  parallel  to  the  x  axis  and 
computed  the  losses  at  different  x  stations.  These  values  are 
plotted  in  Fig.  7.  The  coordinates  of  the  midpoint  between  a 
maximum  and  a  minimum  are  shown  by  the  crosses  in  Fig. 
8,  indicating  the  location  of  an  oblique  shock  at  an  angle  of 
about  41  deg  with  the  x  axis  due  to  the  compression  im¬ 
mediately  after  the  PrandtI-Meyer  expansion. 

A  close  look  at  the  velocity  field  in  Fig.  9  reveals  the 
qualitative  structure  of  the  flow  at  the  shoulder,  which  is  of 
prime  interest.  By  measuring  the  flow  angle,  it  is  found  that 
the  flow  overexpands  as  much  as  -J.I  deg  at  a  point  0.3 
chord  above  and  0.6  chord  downstream  from  the  shoulder. 
A  local  analysis'^  shows  that  the  flow  expands  to  a  max¬ 
imum  velocity 

u|r  .0)-Ar-t-(3/2)»'>-1.8I0 

right  after  the  shoulder  and  compresses  along  the  wall  at  the 
rate 

u-(3/2p^(l-2-»'’3x^’)-t-Ar 

which  is  quite  close  to  the  numerical  prediction.  Because  of 
the  compression,  the  characteristics  coalesce  into  a  weak 
shock  wave  staning  with  a  slope  of  (2/3)*'’  w41  deg,  a  zero 
curvature,  and  a  negative  third  derivative.  The  maximum 
velocity,  linearly  extrapolated  from  downstream  to  the 
shoulder,  assumes  the  value  1.878  on  the  base  grid  and  the 
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values  1.864  and  1.819  after  the  Ttrst  and  second  cycles  of 
reflnemem,  respectively,  approaching  the  theoreiicaj  value  as 
the  solution  on  the  base  grid  is  improved.  It  is  this  oblique 
shock  that  accounts  for  the  balance  of  momentum. 

CondiuioM 

An  efftcient  refinement  procedure  combining  the  ideas  of 
solving  a  modified  difference  equation  and  of  adaptive  mesh 
refinement  has  been  developed  for  refined  computations  of 
steady,  inviscid,  transonic  flows.  The  advantage  of  this  pro¬ 
cedure  is  that  it  allows  local,  as  well  as  global,  improvement 
on  a  base  solution  without  changing  the  base  grid  structure 
and  with  only  modest  increases  in  storage  and  computation 
time. 

In  the  computation  of  the  flow  over  a  wedge  at  AfwO.S, 
the  existence  of  an  oblique  shock  near  the  shoulder  is  con¬ 
firmed  by  comparing  the  pressure  drag  and  the  wave  drag. 

Although  the  present  method  is  applied  to  obtain  solutions 
of  the  small-disturbance  equation  for  a  particular  profile,  it 
can  be  extended,  in  principle,  to  other  complicated  flow 
problems. 
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In  the  context  of  finite  difference  approximation,  the  difference  between  a  differential  equation  and 
the  discrete  (difference)  equation  that  models  it  is  the  truncation  error,  and  it  can  be  easily  shown  that 
if  the  truncation  error  were  known,  solving  the  discrete  equation  would  yield  the  exact  solution  at 
nodal  points.  An  adaptive  procedure  for  improving  the  accuracy  of  a  numerical  solution  on  a  fixed  grid, 
which  we  call  the  base,  through  the  approximation  of  the  truncation  error  by  subdomain  gnd 
refinements  is  introduced.  Regions  where  refinements  are  needed  are  identified  using  an  estimate  of 
the  truncation  error.  Local  solutions  on  grids  tailored  to  each  of  the  regions  are  constructed  and  used 
to  form  approximations  to  the  truncation  error,  which  is  then  'injected'  into  the  base  grid  to  improve 
the  base  solution.  A  one>dimensional  model  of  the  convection>diffusion  equation  is  used  to  demon¬ 
strate  the  basic  ideas  behind  this  method;  two  other  e.xample3  which  imply  extensions  of  this  method  to 
multidimensional  problems  are  also  studied. 


1.  Intrtxluctioa 

Since  the  advent  of  the  computer,  scientists  and  engineers  have  been  using  discrete 
approximations  to  obtain  solutions  to  differential  equations  of  various  complexities.  Most 
problems  can  be  solved  using  a  computer,  given  enough  storage  and  computation  time.  In 
practice,  the  class  of  problems  that  can  be  solved  is  directly  limited  by  the  available  memory 
capacity  and  speed  and,  indirectly,  by  the  efficiency  of  the  method  used  to  calculate  the 
solution.  In  this  paper,  we  introduce  a  grid-refinement  procedure  which  should  substantially 
improve  the  efficiency  of  obtaining  uniformly  accurate  numerical  solutions. 

In  solving  a  differential  equation  that  models  a  physical  problem  of  interest,  e.g.,  the 
Navier-Stokes  equations  in  fluid  mechanics,  the  solution  is  represented  by  values  at  points 
distributed  over  the  space  on  which  the  problem  is  defined.  This  distribution  of  grid  points  is 
usually  at  one’s  discretion  and  is  generally  related  to  the  geometry  and  some  characteristic 
features  of  the  problem.  It  directly  affects  the  accuracy  of  the  solution.  The  optional 
distribution  of  grid  points  for  various  problems  has  been  attempted  by  many  researchers.  A 
review  by  Thacker  [1]  surveyed  the  various  methods  currently  used  to  generate  grids  that  are 
geometry  related.  These  methods  are  concerned  with  the  regularity  and  smoothness  of  the 
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distribution  of  the  grid  generated  and  the  satisfaction  of  certain  predetermined  criteria  that  are 
believed  to  be  essential  for  accurate  solutions.  However,  there  is  no  guarantee  that  the 
solution  obtained  with  such  grids  will  meet  cenain  accuracy  requirements,  nor  that  the 
number  of  grid  points  will  be  optimum  for  a  given  accuracy;  such  goals  can  only  be  realized  by 
grids  that  are,  at  least  to  some  degree,  solution  dependent. 

As  pointed  out  in  a  recent  paper  by  Benek,  Steger  and  Dougherty  [2],  solution-dependent 
grid-generation  methods  can  be  classified  into  three  categories:  grid  patching,  grid  embedding, 
and  grid  adapting.  The  applications  of  these  methods  are  problem  oriented.  Grid  patching  and 
grid  embedding  are  based  on  the  idea  of  dividing  the  problem  domain  into  subdomains  which 
are  easier  to  handle  or  over  which  the  solutions  can  be  obtained  more  efficiently  than  over  the 
full  problem.  A  typical  example  is  the  aircraft  design  and  analysis  problem  [3,4].  Grid-to-grid 
communication  is  handled  by  interpolations,  and  care  must  be  taken  at  intergrid  boundaries  to 
ensure  certain  compatibility  conditions  are  satisfied,  e.g.,  conservation  of  fluxes  [5].  Usually, 
this  type  of  grid  generation  is  only  problem  dependent,  designed  according  to  the  user’s  view 
of  the  problem,  but  can  be  made  flexible  enough  to  move  around  the  problem  domain,  change 
size  and  density,  and  retire  if  no  longer  needed  [6].  These  procedures  require  extensive 
bookkeeping  for  both  the  grid  system  and  the  solution. 

There  are  two  grid-adapting  strategies.  One  strategy  is  to  redistribute  a  fixed  number  of  grid 
points  according  to  some  criteria  such  that  the  overall  error  measured  by  some  means  is 
reduced  (e.g.,  [7, 8]).  The  other  is  to  increase  the  number  of  grid  points  near  regions  where  a 
measure  of  error  indicates  refinement  is  needed  (e.g.,  [9, 10]).  These  methods  are  most 
effective  for  problems  with  local  singularities  or  layers  of  rapid  changes,  e  g.,  shock  and 
boundary  layers.  However,  a  successful  application  of  these  methods  depends  on  the  criteria 
used  for  grid  redistribution  or  addition,  which  may  be  very  sensitive  to  the  behavior  of  the 
solution,  and  repeated  computations  are  needed  before  a  satisfactory  result  is  achieved.  In  the 
case  of  grid  addition,  redistribution  of  grid  points  may  also  be  needed  to  avoid  the  error  due 
to  abrupt  changes  of  grid-point  distribution.  For  multidimensional  problems,  the  extra  effort 
needed  to  manage  the  complexities  associated  with  adaptive  gridding  could  upset  the  overall 
effectiveness  of  these  methods.  Advanced  computer  science  techniques  like  pattern  recogni¬ 
tion,  artificial  intelligence,  and  data  base  management  will  definitely  play  a  significant  role  in 
reducing  the  total  cost  of  obtaining  a  numerical  solution. 

A  natural  question  arises  in  refining  numerical  solutions:  is  it  necessary  to  increase  the 
number  of  grid  points  or  redistribute  them  in  order  to  improve  the  accuracy  of  a  numerical 
solution?  The  answer  is  not  as  straightforward  as  one  might  think.  Customarily,  the  error  of  a 
numerical  solution  is  related  to  the  truncation  error  or  the  residual  due  to  discretization. 
When  a  measure  of  the  truncation  error  is  larger  than  a  preset  value  in  certain  regions,  the 
addition  of  sufficient  grid  points  in  those  regions  reduces  the  truncation  error  and,  hence, 
improves  the  accuracy  of  the  solution,  a  consistency  requirement.  However,  it  will  be  shown 
in  Section  2  that  an  improved  solution  is  possible  without  changing  the  base  grid  structure  by 
using  an  approximation  of  the  truncation  error  obtained  from  a  local  refinement.  These  two 
approaches,  while  similar  in  nature,  are  different  in  methodology. 

An  adaptive  refinement  procedure  for  solving  differential  equations  is  introduced  in  Section 
3.  Applications  of  this  procedure  to  typical  problems  in  fluid  mechanics  are  demonstrated  in 
Section  4. 
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2.  Truncation  error  and  grid  refinement 

Most  numerical  analysis  books  state  that  a  differential  operation  L  operating  on  a  function 
(p  is  related  to  a  difference  operator  (with  the  subscript  h  denoting  the  grid  size),  operating 
on  (p  by  the  truncation  error  (abbreviated  TE  hereafter),  i.e., 

Ltp  *L^(p  +TE(<p,  h) .  (1) 

This  result  of  a  Taylor-series  expansion  is  the  basis  for  all  finite  difference  techniques. 
Ordinarily,  the  direct  system, 

(2) 

is  solved  for  tp^,  which  satisfies,  consequently,  the  equation 

L<p^=TE(<p^,h) .  (3) 

As  a  consistency  requirement,  the  TE  vanishes  as  the  step  size  h  approaches  zero,  leading  to 
the  limiting  solution  <p, 

lim  , 

h-*0 

that  satisfies  the  differential  equation 

Lip^O.  (4) 

Hence,  it  is  assumed  that  if  the  TE  is  uniformly  small,  solving  the  discrete  system  (2)  will 
lead  to  a  good  approximate  solution  of  <p.  However,  in  many  problems  the  TE  is  a  rapidly 
varying  function  of  its  arguments.  The  idea  of  conventional  grid-adapting  techniques  is  to  look 
for  or  modify  the  distribution  of  grid  points  according  to  some  preset  criteria  which  will  render 
uniformity  of  truncation  error  across  the  solution  domain.  Unfortunately,  there  is  no  simple 
way  to  generate  a  grid  that  minimizes  the  truncation  error  for  a  given  problem  and,  in  many 
cases,  the  process  of  finding  the  optimal  grid  is  more  complicated  and  time  consuming  than 
computing  the  solution  itself. 

We  must  remember  that  it  is  (4)  that  one  wants  to  solve,  not  (3),  which  is  equivalent  to 
solving  the  discrete  equation  (2).  The  discrete  equation  that  ought  to  be  solved  corresponding 
to  solving  (4)  is  implied  by  (1),  i.e., 

+TE((P, /|)»0.  (5) 

Here,  we  have  deliberately  denoted  the  argument  function  9  of  the  TE  with  a  tilde,  which  can 
be  different  from  the  solution  <p^.  Notice  that  if  the  exact  solution  were  available,  it  would 
satisfy  (5)  exactly,  with  (p  *  **  ip.  This  implies  that  the  TE  can  be  computed  exactly  by 

applying  the  operator  L^,  to  the  solution,  e.g.. 
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TE((p,  -L*<p  ,  (6) 

and  that  solving  the  system 

h)  =  L^<p 

will  yield  the  exact  solution  at  nodal  points.  Hence,  it  is  clear  that  if  the  goal  is  to  improve  a 
numerical  solution,  the  basic  grid  structure  need  not  be  changed,  only  improved  values  of  TE 
at  grid  points  need  to  be  provided.  To  emphasize  this  point,  the  base  grid  in  the  examples 
considered  here  is  never  changed.  The  strategy  one  would  use  in  making  tradeoffs  between 
the  base  grid  and  approximating  the  TE  is  not  discussed  in  this  paper. 

Analytically,  the  TE  consists  of  all  higher  derivatives  of  the  function  being  expanded  in  a 
Taylor  series,  if  they  exist.  If  more  neighboring  values  of  a  function  are  known,  higher 
derivatives  can  be  computed,  and  hence  the  TE  can  be  better  approximated.  A  TE  sequence 
may  be  defined  as  follows: 

TE/»/,v*TE(^*,,v.^)* (7) 

where  the  subscript  hJN  refers  to  values  based  on  a  grid  of  size  h/N  (e.g..  subdividing  the  base 
grid  of  size  h,  N  times).  For  simplicity,  we  may  assume  satisfies  the  equation 

^hlN^hlS  •  (8) 

With  these  definitions  and  the  substitution  of  (7)  into  (6)  and  setting  ^  v»  is  easy  to  see 
that  <Pf,  =  <p^,^  is  a  solution  of  (5)  at  coinciding  nodal  points  (or  through  the  use  of  an 
interpolation  function). 

All  we  have  shown  so  far  is  that  it  is  possible  to  obtain  a  refined  numerical  solution 
satisfying  (5)  without  changing  the  base  grid,  provided  the  TE  is  known  to  the  same  order  of 
accuracy  as  the  solution.  The  inclusion  of  TE  into  the  difference  equation  had  been  suggested 
before  [11-14].  Pereyra  [11]  proposed  a  deferred  correction  procedure  to  compute  the  TE 
progressively.  Warming  and  Hyett  [12]  and  Klopfer  and  McRae  [13]  implemented  it  in  forms 
of  a  modified  equation.  For  some  simple  linear  differential  equations,  even  analytical 
expressions  of  the  truncation  error  in  terms  of  lower  derivatives  of  the  unknown  function  have 
been  used.  However,  the  complexity  involved  in  the  derivation  of  such  terms  and,  in  some 
cases,  the  numerical  instabilities  caused  by  the  presence  of  certain  terms  has  discouraged  the 
more  popular  use  of  such  schemes.  In  the  following  sections  we  will  introduce  an  adaptive 
procedure  similar  to  that  suggested  by  Brandt  [14]  for  solution  refinements  based  on  the  ideas 
discussed  in  this  section  and  show  that  the  injection  of  the  TE  is  a  simple  but  effective  means 
of  improving  the  accuracy  of  a  numerical  solution. 


3.  A  refinement  procedure 

Boundary  and  internal  layers  are  common  structures  in  nonlinear  mechanics.  Analytically, 
such  layers  are  commonly  found  through  singular  penurbations,  whereby  solutions  in  different 
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regions  are  described  by  simplifed  versions  of  the  governing  equation  obtained  using  a  scaling 
pertinent  to  the  local  features.  The  grid  rehnement  procedure  used  here  is,  in  many  ways, 
analogous  to  the  singular  perturbation  method.  Clustering  grid  points  over  areas  where  rapid 
changes  occur  is  effectively  scaling  the  local  solutions.  A  discrete  operator  can  be  regarded 
as  a  simplified  operator  of  the  full  equation  and  is  valid  everywhere  e.xcept  in  regions  where 
the  omitted  higher-order  terms,  TE(^,  h),  become  dominant. 

The  refinement  procedure  we  have  chosen  uses  the  TE  as  an  indicator  and  as  a  means  for 
improving  the  solution.  To  begin,  a  base  grid  of  size  h,  not  necessarily  uniform,  is  chosen  for 
the  base  solution  <p^.  The  TE  is  initially  assumed  to  be  zero.  Equation  (5)  is  then  solved  for 
to  a  preset  accuracy  e.  Regarding  as  a.  refined  solution  for  the  grid  2h,  which  is  formed  by 
omitting  every  other  grid  point  of  the  base  grid,  the  truncation  error  is  then  estimated  using 
the  formula  for  every  point  of  the  grid  h  except  for  the  boundary  points  and  the  ones 
next  to  them.  Tlie  regions  where  the  estimated  TE,  is  larger  than  a  preset  value  e  (or 

other  appropriate  criterion)  are  identified;  for  this,  a  pattern  recognition  algorithm  like  that 
described  in  [6]  wilt  be  very  useful.  Outside  the  region  where  the  truncation  error  is  injected. 
TE  is  set  to  zero.  Buffer  zones  are  then  introduced  with  TE  »  0  but  with  the  mesh  refined  to 
achieve  smooth  transition  of  solutions.  Information  on  the  size  of  the  region  with  TE 
injections,  solution  boundary  values,  and  parameters  like  the  refinement  factor  N  of  these 
regions  are  then  transmitted  to  the  finer  grid  solver,  I.a/,v«’a/.v  obtain  a  refined  local 

solution,  v*/,v  These  local  refined  solutions  are  then  used  to  form  the  approximated  TE. 
for  the  next  cycle  of  refinement  until 

®  fo^  M  >  N . 


AOAPTive  CRiO  RCFINEMCNTS 


Fig.  1.  Procedure  for  adaptive  grid  refinement  with  truncation  error  injection. 
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At  this  stage,  further  refinements  will  have  no  effect  on  the  solution  ip,,  which  then  satisfies  (S) 
to  the  preset  accuracy  e.  This  procedure  is  shown  schematically  in  Fig.  1. 

Note  that  each  refinement  may  have  its  own  local  solutions  injecting  their  TE  either  to  their 
parent  solutions  or  to  the  base  solution,  its  own  governing  equations,  and  its  own  grid 
distribution.  No  attempt  is  made  here  to  optimize  the  strategy  for  grid  nesting.  In  the  case  of 
using  a  local  stretched  grid,  interpolation  is  needed  to  obtain  the  TE.  A  major  advantage  of 
truncation  error  injection  is  that,  unlike  other  grid-refinement  methods,  there  is  no  particular 
need  for  storing  the  finer-grid  solutions.  Once  the  TE  is  obtained,  the  approximated  values  of 
the  TE  indicate  how  good  the  local  solution  is  and  if  needed  refinements  can  be  obtained  with 
minor  effort. 


4.  Numerical  examples 

The  first  example  chosen  to  demonstrate  the  adaptive  procedure  was  the  classical  example 
in  boundary  layer  theory  modeled  by  the  one-dimensional  convection-diffusion  equation: 

<0,  *  ^  ,  with  (p(0)  “0  and  (p(l) «  I . 

This  example  has  been  used  by  many  authors  to  demonstrate  special  properties  of  their 
methods. 

In  particular,  we  chose  central  differences  for  both  the  first  and  second  derivatives.  It  is  well 

<9  h»l/8 
- EXACT  ip 


X 


REYNOLDS  NUMBER  •  100 
CENTRAL  DIFFERENCE  SCHEME 
SOLUTIONSOF  hX/RE-oV®**  *0 

Fig.  2.  Solution  without  TB  injection. 
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known  that  for  a  grid  Reynolds  number  Re  >  2,  nonphysical  spatial  oscillations  due  to  large 
dispersive  truncation  error  occtir.  Figure  2  shows  the  erratic  solution  for  Re  » 100  and  step 
size  /t  *  i .  It  bears  no  resemblance  to  the  exact  solution,  which  remains  flat  until  near  the 
right  boundary  and  shoots  up  exponentially  to  the  boundary  value  ^(1)  1. 

After  a  TE  estimation  (dashed  line  in  Fig.  3(a))  which  indicates  that  refinement  is  needed 
everywhere,  the  chosen  subdomain  (same  as  the  base  grid)  is  divided  so  it  has  a  local  step  size 
of  ih.  A  refined  soution  (p^,2  is  then  obtained,  solving  (8).  and  the  corresponding  is 
computed  and  injected  into  (S).  The  solution  obtained  and  absolute  error  are  depicted  in  Figs. 
3(b),  (c).  The  effects  of  the  truncation  error  are  evident  by  the  disappearance  of  the 
nonphysical  oscillations  in  the  recomputed  solution  (p^.  Progressive  improvements  of  the 
solution  after  the  second  and  third  cycle  of  refinement  can  be  seen  in  Figs.  4  and  3.  Because  of 
the  preset  tolerance  (e  »  10~^),  the  refined  region  was  the  whole  region  for  the  above  cases. 
Figure  6  shows  a  converged  solution  accurate  within  10~^  in  absolute  error  after  the  fourth 
refinement  cycle.  From  the  estimated  TE,  the  algorithm  decided  that  it  was  sufficient  to  refine 
the  grid  from  x  >  0.373  to  x  »  1.  The  total  number  of  grid  points  used  in  the  last  case  was  SI. 
only  a  few  points  more  than  the  previous  .eflnement.  It  is  interesting  to  note  that  in  the 
converged  solution  the  error  is  uniformly  distributed  across  the  solution  domain.  Further 
refinements  showed  no  further  improvement  in  the  solution,  indicating  that  the  last  local 
refinement  was  enough  to  satisfy  the  preset  error  criterion.  Since  this  is  a  linear  problem,  a 
directed  computation  using  217  uniformly  distributed  points,  the  same  total  number  of  points 
used  in  this  procedure,  shows  a  total  mean*squared  error  of  0.023;  an  improvement  of  10'^  is 
achieved  using  this  refinement  strategy  with  the  same  amount  of  computational  resources. 

It  is  of  interest  to  note  that  the  choice  of  a  scheme  is  immaterial  to  the  effectiveness  of  this 
solution  refinement  method.  We  have  applied  it  to  the  same  problem  using  an  upwind  scheme. 


- CSTIM4TE0  TE 

- APPBOXIMATCO  TE 


O  xith  TE  PPOM  «  U6PTS 

n/c 


RCPINC  LOCALLY  IP  TE*  .0001 
PROM  POINT  I  TO  POINT  9 
ESTIMATED  TRUNCATION  ERROR 


CONVEROEO  IF  CMANCES  0P^<  0001 
ERROR  OP  SOLUTIONS 


REYNOLDS  NUMBER  •  lOO 
CENTRAL  DIFFERENCE  SCHEME 
SOLUTIONS  OF  liX/RE-0**/0X*-0^X»O 


Fig.  3.  Solution  with  TE  injection  after  first  cycle. 
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- CSTIMATCO  TE 

- APPROXIMATED  TE 


«  WITH  TE  FROM  132  PTS) 


REFINE  locally  iF  TE  >  .0001  CONVEROCO  IF  CHANGES  OF  ^  <  OCOI  REYNOLDS  NUMdER  •  100 

FROM  POINT  I  TO  POINT  9  CENTRAL  DIFFERENCE  SCHEME 

ESTIMATED  TRUNCATION  ERROR  ERROR  CF  SOLUTICNS  SOLUTIONS  OF  h/RE-  9*^/OX*‘-0(^OX>0 

Fig.  4.  Solution  with  TE  injection  after  the  second  cycle. 


The  results  are  even  more  impressive.  A  solution  accurate  to  10 is  achieved  with  two  levels 
of  local  refinement  and  a  total  effort  equivalent  to  that  of  using  29  points.  However,  the 
examples  using  central  differences  are  more  illustrative. 

The  second  example  was  the  reflections  of  a  two^imensional  oblique  shock  from  a  wall.  We 
solve  the  Euler’s  equation  using  the  MacCormack's  scheme  because  of  its  simplicity  and 
popularity.  Figure  7  shows  a  display  of  the  exact  pressure  distribution  for  an  incident  shock,  a 


- estimated  TE  TE  FROM  ^,,3  (64  PTS) 

—  APPROXIMATED  TE 


TE  E 


0.00  0.23  0.30  0.T3  1.00  0.00  023  000  0.73  1.00  0  00  0  23  0.30  073  100 

XXX 

So  ot  Sc 


REFINE  LOCALLY  IF  TE  >.0001  CONVERGED  IF  CHANGES  OF  ^  <.0001  REYNOLDS  NUMBER  •  100 

FROM  POINT  I  TO  POINT  9  CENTRAL  DIFFERENCE  SCHEME 

ESTIMATED  TRUNCATION  EPflOR  E’ROR  OF  SOLUTION”.  SOLUTIONS  OF  h/RE- 0*4/0X*-0^/0X.C 

Fig.  5.  Solution  with  TE  injection  after  the  third  cycle. 
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Fig.  7(a).  Pressure  contours  of  an  oblique  shock  reflection,  (b)  Exact  truncation  error  computed.  TE  values  have 
an  increment  value  of  0.001,  starting  at  a  base  value  of  0.01. 
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pressure  discontinuity  of  35*  at  free  stream  Mach  number  =  2.013  over  a  11  x  31  grid  of 
step  size  h  »  0.05.  Fifteen  contour  levels  are  plotted  in  increments  of  0.5  beginning  at  0.73. 
The  zigzags  are  a  result  of  representing  on  a  grid  the  abrupt  pressure  jump  from  one  nodal 
point  to  another,  depending  on  whether  the  point  is  before  or  after  the  shock,  by  interpolated 
contours  for  graphic  display.  For  reference,  the  exact  truncation  error  computed  from  the 
exact  solution  (Fig.  7(a))  is  shown  in  Fig.  7(b)  in  increments  of  0.005  starting  at  0.001.  It  is 
clear  from  Fig.  7  that  local  refinements  with  rotated  grids  aligned  with  the  shock  would  be 
most  effective.  However,  our  intention  in  this  example  was  to  show  that  even  for  a  nonlinear 
operator,  the  exact  solution  can  be  recovered  if  the  TE  is  given  or  computed  locally.  Figure  8 
shows  the  approximated  truncation  error  from  a  refined  solution  with  a  refinement  factor  of 
four.  A  comparison  of  the  pressure  contours  before  and  after  TE  injection  shows  marked 
improvement.  Except  for  the  small  wake  behind  the  shocks,  which  is  due  to  the  numerical 
scheme  selected,  the  refined  solution  is  very  close  to  the  exact  one. 

The  third  problem  was  chosen  to  study  this  procedure  with  rotated  grids  aligned  with  the 
discontinuities.  We  solved  the  linear,  two-dimensional  convection-diffusion  equation  for  a 
nominal  quantity  7,  modeling  two  adjacent  fluids  of  initially  different  temperatures  moving  at 
the  same  speed.  Upwind  differencing  was  used  for  the  convective  terms.  It  is  well  known  that 
the  artificial  cross-wind  diffusion  due  to  upwind  differencing  is  a  major  source  of  error.  It 
causes  excessive  spreading  of  the  discontinuity.  Figure  9  shows  the  temperature  contours 


TRUNCATION  ERROR.  RF«4 


TRUNCATION  ERROR  CONTOURS  COMPUTED 
FROM  A  REFINED  SOLUTION  WITH  A 
REFINEMENT  FACTOR  OF  4. 


SEFORE  l-„4h*0 


Fig.  8.  Pressure  contours  before  and  after  truncation  error  injection,  W  -  4. 
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Fig.  10.  Numerical  solution  of  a  temperature  shear  layer,  grid  120  x  30,  no  TE  injection. 

solved  on  a  60  x  40  base  grid  of  step  size  h  -  0.25.  A  slightly  improved  solution  on  an 
unrotated  120  x  80  grid  is  shown  in  Fig.  10  for  comparison.  This  does  not  compare  well  with 
the  exact  solution  depicted  in  Fig.  11.  Given  suitable  pattern  recognition  schemes,  it  would  be 
natural  to  introduce  a  rotated  grid  parallel  to  the  flow  direction  over  a  small  region 
surrounding  the  discontinuity  with  boundary  conditions  extracted  from  the  solution  in  Fig.  9. 
Here,  we  have  done  this  manually  with  a  20  x  40  grid  aligned  with  the  flow.  Due  to  grid 
rotation,  the  cross-wind  artificial  diffusion  is  minimized,  resulting  in  a  sharp  temperature 
gradient  very  close  to  the  exact  solution.  The  isotherms  that  appear  near  the  upper  and  lower 
boundary  of  the  refined  local  solution  on  the  subgrid  (Fig.  12)  are  an  effect  of  the  incorrect 
boundary  conditions  extracted  from  the  base  grid  solution;  these  can  be  avoided  simply  by 
taking  a  larger  subgrid.  However,  with  the  injected  TE,  the  improved  base  solution  provides  a 
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sharp  gradient  without  such  isotherms  (Fig.  13),  which  shows  that  the  base  grid  solution  is 
readjusted  smoothly  through  the  artificial  boundaries. 

The  examples  given  here  demonstrate  various  aspects  of  the  TE  injection  method  and 
show,  in  particular,  that  very  accurate  solutions  can  be  obtained  on  relatively  coarse  meshes. 


5.  Conclusions 

It  is  well  known  that  if  the  truncation  error  in  a  discrete  approximation  to  a  differential 
equation  were  known  exactly,  the  values  of  the  exact  solution  at  the  discrete  points  could  be 
determined.  This  fact,  along  with  adaptive  mesh  refinements  to  determine  the  truncation 
error,  is  used  to  produce  highly  accurate  solutions  to  model  problems  on  relatively  coarse 
grids  with  substantial  savings  in  computer  time  and  use  of  computer  memory.  Improvements 
in  the  local  as  well  as  global  accuracy  of  a  solution  on  a  fixed  grid  are  found  by  refining  the 
grid  to  estimate  the  truncation  error  and  by  injecting  this  truncation  error  back  into  the 
solution  of  the  discrete  equation  on  the  unrefined  grid.  This  use  of  solution  adaptive  grid 
systems  reduces  the  dependence  of  the  solution  on  the  choice  of  the  grid  and,  hence,  the  effort 
of  grid  generation. 
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